function [XBoot, RmBoot, RfBoot] = BootGenDataAR1(X, Rm, Rf, NBoot)
% DGP is assumed to be AR(1) for X and iid for Rm and Rf

[T, K] = size(X);
RmBoot = NaN(T, NBoot); RfBoot = NaN(T, NBoot); XBoot = NaN(T, NBoot, K);

% If X is a matrix, then it can happen that it is not a full panel.
% Meaning some values X(t,k) are missing. The code handles this.
% And generates the bootstrapped sample that is also NOT full panel.

%% Generate in-sample estimates
muRm = mean(Rm); eRm = Rm - muRm;
muRf = mean(Rf); eRf = Rf - muRf;
muX = mean(X,'omitnan'); bX = NaN(2,K); eX = NaN(T,K);
for k = 1:K
    Xk = X(:,k);
    y = Xk(2:end); x = [ones(T-1,1) Xk(1:end-1)];
    Tidxk = isfinite(y) & all(isfinite(x),2);
    x = x(Tidxk,:); y = y(Tidxk);
    b = x\y;
    b(2) = b(2) + (1+3*b(2))/(sum(Tidxk) - 1); %Correct for small-sample bias of AR(1) coefficients
    if b(2)>=0.999; b(2) = 0.999; end
    b(1) = muX(k)*(1 - b(2));
    bX(:,k) = b;
    e = NaN(T-1,1); e(Tidxk) = y - x*b;
    eX(2:end,k) = e;
end

%% Generate bootstrapped data
Textra = ceil(T/10); % these are the 10% extra observations to generate
% Now generate NBoot*(T+Textra) random numbers between 2 and T
TIndex = randi([2 T],NBoot*(T+Textra),1); TIndex = reshape(TIndex,T+Textra,NBoot);
for b = 1:NBoot
    Tb = TIndex(:,b);
    eRmb = eRm(Tb); eRfb = eRf(Tb); eXb = eX(Tb,:);
    Rmb = muRm + eRmb; RmBoot(:,b) = Rmb(Textra+1:end);
    Rfb = muRf + eRfb; RfBoot(:,b) = Rfb(Textra+1:end);
    for k = 1:K
        Xb = muX(k)*ones(T+Textra,1); % set 1st observation to mean
        tprev = 1;
        for t = 2:T+Textra
            eXb_tk = eXb(t,k);
            if isfinite(eXb_tk)  % need to do this to handle missing eXb
                Xb(t) = bX(1,k) + bX(2,k)*Xb(tprev) + eXb_tk;
                tprev = t;
            else
                Xb(t) = NaN;
            end
        end
        XBoot(:,b,k) = Xb(Textra+1:end);
    end
end
return